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The measurement process uncertainty is propagated through the use of a calibration curve. The magnitude 
and direction of this uncertainty depends on the choice of the controllable variable in producing the calibration 
curve; in other words, the design of the calibration experiment. In this paper this design is discussed in the con- 
text of Scheffe's approach to the uncertainties of a calibration curve and in particular for the case in which the 
calibration curve is a linear spline. A class of appropriate designs is given, which depend on the location of the 
knots and the slopes of the segments. One of these designs is quickly calculable and can be found without a com- 
puter. Based on these results, a design approach is suggested for the case in which the knots are not known 
exactly. 
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1. Introduction 

This paper considers some design aspects of the calibration problem within the uncertainty formulation 
given by Scheffe [1973]. A mathematical and statistical relationship is assumed to exist between two quanti- 
ties U and v where v is generally more expensive or difficult to measure than U. For a given value oft?, obser- 
vations on U are assumed to be random with a mean value 

m(v) = m(v,p) = iPtgiv) (1.1) 

and variance o 2 independent of v. Of course the functions g if t = 0,1, . . . ,A;-F 1 form the basis for the regres- 
sion function m, and the scalars /?,, i=0, . . . ,i+ 1 and o are unknown. For given values [/,*, C/ 2 *, . • . , one is 
interested in finding the corresponding v x *, v 2 *> .... To do this the system is "calibrated." That is, exact 
values v iy v* . . . >v n are chosen for v, and corresponding readings U u U^ . . . ,U n are taken. The regression 
coefficients ($ are then estimated by the least squares estimates, ft. For given values U* - u* the corre- 
sponding estimates oft;,*, v h are found by solving u* = m(v iy p). 

Scheffe [1973] shows how a "calibration chart" can be constructed which results in a band around the 
curve u = m(v,p). (Throughout we shall view u and v as plotted in the vertical and horizontal directions 
respectively, and thus differ from Scheffe by using the conventional axes.) These bands are used to construct 
an interval estimate I(u) for the unknown value of v. The bands are constructed for a given a and d so that 
"for every possible sequence of constants v* the probability is at least 1—6 that the proportion of intervals 
containing the corresponding v t * is in the long run at least 1— a." The reader is referred to Scheffe [1973] 
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and Scheffe, Rosenblatt, and Spiegelman [1980] for further discussion and interpretation of the bands. 
Alternative approaches to the uncertainty of the calibration problem can also be found in Lieberman, 
Miller, and Hamilton (1967) and the National Bureau of Standards Special Publication 300 on Precision 
Measurement and Calibration, Ku (1967). 

In order to prevent the diversion (theft) of nuclear material it is important to be able to determine, within 
the context of statistical uncertainties, the amount of radioactive material in the tank when a pressure 
measurement is taken. It is clear that the larger the uncertainty limits on the volume estimates, the more 
likely it is that the diversion of some material will go undetected. For this reason designs which result in 
large statistical uncertainties for volume estimates are to be avoided. Therefore, the designs for calibration 
experiments given here seek to minimize the maximum uncertainty limits that result from the use of a cali- 
bration curve. 

In order to obtain the shape of the calibration curve, m, the basic relationships of physics between pres- 
sure and volume are given. For a homogeneous liquid the pressure reading is proportional to the height of 
the liquid in the tank, h. Since volume — v — J o A(x)dx, where A(x) is the cross sectional area at height x, and 
pressure = P = gdh, where g is the acceleration due to gravity, and d is the density of the liquid, the 

P/gd 

volume-pressure relationship may be written as v = J A(x)dx. Thus the pressure-volume relationship of a 
constant density liquid is a straight line if the cross sectional area^(x) is constant 

However often there are significant obstructions in the tank. In these cases the cross sectional area is not 
constant and the calibration relationship is not well-approximated by a straight line. These obstructions are 
of the nature of cooling coils, mechanisms to agitate the liquid or supporting metal to strengthen the tank. 
An abrupt obstruction would cause a sharp change in slope of the pressure volume relationship. When all 
obstructions are abrupt, relative to the size of the tank, it can be seen that a linear spline (given below) is 
appropriate. (In some cases obstructions would change the cross sectional area in a significantly more 
gradual manner, in these cases nonlinear splines might be appropriate.) 

In this case the pressure volume relationship can be usefully represented as a linear spline 



m(v)=a+ bv +1(^)^-1) 



(1.2) 



where z + = max {0,z}. The quantities |i <| 2 < . . . <!* are called f "knots" and these will occur wherever the 
cross-sectional area changes. In eq (1.2) the function is a + bv for v below f t . The curve is continuous and 
changes to slope b + b x at the point % u etc. 

When calculating the estimates of the calibration curve from the data it is often convenient to work with a 
more complicated but also more numerically stable basis than that used in (1.2). We use | and |* +1 for the 
smallest and largest volume readings respectively. Then define 



NM = 



N,{v) 



l.-lo 


1,-1.-. 





lo<f<l, 
elsewhere 



for i = 1,2,. . ,,k 



N k .M = 







!*< v<!* + i 



»>!*♦! 



These are called 5-splines. They are simply triangular type functions over successive pairs of intervals. With 
this basis the value /J, in 



296 



m(v) 



Ipmv) 



(1.3) 



is given by the value of m(v) at |,, i.e. /?, = tfi(| t ), i = 0, 1, . . . ,A+ 1. An excellent discussion of splines is 
given in de Boor [1978]. 

The Scheffe calibration bands around the function (1.3) produce various width inverse intervals I(u) for 
the corresponding unknown value v. It is clear that the wider the band is in the vertical u direction the wider 
are the horizontal intervals I(u). In addition, however, low slope values of the curve m(v) in (1.3) will produce 
much larger intervals I(u) than large slopes will. 

For the reason explained earlier we would like to have the maximum (over the range of u values) of the 
lengths of the intervals I(u) as small as possible. This will also be useful when one wishes to make a single 
quantitive statement about the overall accuracy of the calibration. It is possible to keep the maximum of I(u) 
small if higher precision can be obtained in estimating the true curve v(v) in regions where it has lower 
slope. This precision is measured through the variance or standard deviation of our estimate of the response 
curve m{v). More observations are needed where the knots occur and where the slopes are low. The design 
problem is to see what quantative statements can be made about where the values v u v* . . . ,v n should be 
chosen to obtain corresponding readings U u U* . . . ,U n in the calibration experiment 

Typically the first few values of v are taken nearly evenly spaced so that a judgment as to whether the 
calibration experiment is in control or not can be made. If the calibration experiment is under control then 
the experimenter can choose the rest of his f V values to minimize, to the extent possible, the uncertainty 
limits that result from using the calibration curve. 

In section 2 some consideration is given to the selection of v u v 2 , . . . ,v n , and a procedure is proposed. An 
illustrative example is described in section 3. Some discussion is given in section 4. 

2. Choosing the volume values 

As mentioned in the previous section it is desirable to keep the maximum width of the intervals I(u) at a 
minimum. In order to do this we consider a general curve (see fig. 1) of the form 



m = m(v,p) = I pMv) 



(2.1) 



where the /3, are all positive and increasing. (The range of v is from | to £ k+l where generally | > 0. At the 
bottom of the tank the volume and pressure should both be zero. However tanks do not have flat bottoms but 
rather are bowed so that the validity of (2.1) is questionable near the bottom. We take | > 0.) 
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Figure 1. 
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In order to minimize max I(u) we draw a band of constant horizontal width 2d about the curve. The 

u 

Scheffe bands are produced by considering 

mfofl)±d[c 1 + c 2 S(v)l (2.2) 

The notation will be explained carefully below. Our intent is to choose d so that the constant width curve 
just contains (2.2). The values v u . . . ,v n will then be chosen so as to minimize d, these values enter mainly 
through the quantity S(v). 

To explain the notation, the quantity 6 is an estimate of the standard deviation o in our pressure readings. 
We will assume for simplicity that o is known and therefore take 6 = o. Some discussion in section 4 will be 
given to this matter. The quantities c x and c 2 are constants depending on the values of a and d mentioned in 
the introduction. The function S(v) is, except for a factor of o, the standard deviation of the estimate of 
m(v y p) for a fixed value v of the volume. If we let 

S*(v) = — W(v) M- 1 W N(v) (2.3) 

where the matrix Miy) is given by M(/i) = / N(v) N'(v) a\i (v). The measure \x is called the design measure [see 
V. V. Federov(1972)] and simply has mass \ln at each observation point v it i = 1,2, . . . ,n(fi = number of 
observations). The observations on U are assumed to be uncorrelated and in general some of the v t values 
could be equal. The elements of M = ilfjfyi) are simply 

m ab = iN a (v t ) N b (v t ), < a, b < k + 1 (2.4) 

The goal is to consider the maximum of the horizontal widths of I(u) of the Scheffe bands and minimize 
this maximum with respect to the values v x , v* . . . ,v H . This seems to be extremely difficult We shall proceed 
by trying to minimize d. 

To hold down the maximum width of I(u) we require that the Scheffe bands (2.2) lie inside the bands 
shown in figure 1. Since the upper and lower bands in figure 1 are m(v+d) and m(v—d) respectively we thus 
require that, for all v, 

m(v + d)- m(v) > o[c x + c 2 S(v)] (2.5) 

m(v) - m(v-d) > o[c x + c 2 S(v)] (2.6) 

It will be shown at the end of this section that o[c x + c 2 S(v)] is convex in v on each segment [|„ £ l+1 ], i = 0, 
1, . . . ,A. The right hand side of (2.5) then consists of convex segments, while the left hand side consists of 
linear segments. Equations (2.5) and (2.6) will then hold provided they hold at the knots of the left hand side. 
Thus we require (2.5) to hold for 

| t , i = 0, 1, . . . ,i+ 1 and |-i i = 1,2,...,* (2.7) 

and (2.6) should hold at the points. 

|„j= 0, l,...,fc+l and& + 4i= 1,2,..., A; (2.8) 

Only points in the interval (| , 4*+i) are considered. We assume that 

€m-&>4 * = 0,l,...,i (2.9) 

then (2.5) and (2.6) will hold if 
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o[c x + c 2 5(| )] < K 

o[c, + c 2 S&)] < k and h t . x i = 1, . . . ,4 (2.10) 

where /i, = ra(|, + d) — m(l<) = m(|^i) — m(4 1+1 — c/). Let 5, denote the slope of ra(i>) on (|„ | l+1 ), i = 0, 
1, . . . ,4. Then since h t = ds { , the requirements (2.10) become 

o[ Cl + c 2 S(l)] < dy< i = 0, 1,.. . ,4+1 (2.11) 

where y = So, y« — min{j to «<_i}, i = 1,... ,4 andy fc+1 = s k . 

Now the design measure or the choice of Vu v* . . . ,v„ enters these equations in S(v), The general problem 
is still to choose the design so that the value of d can be made as small as possible. 

To reduce S(v), observations should be chosen at the known points | , li» . • • ,4*+i- (The case of estimate 
knots |, is considered in section 3.) If we note the dependence of S(v) on \x by S(v,\a) then it is known [see 
W. J. Studden and D. J. Arman(1969)] that for a fixed set of knots and any/i there is another design^,, con- 
centrating on !„ such that S[v,\jl x ) ^ S(v^ ) for all v. 

Suppose, then, that we had observations only at the endpoints and knots. The matrix M( ) in (2.3) and (2.4) 
can be seen to be a diagonal matrix with diagonal elements po, p u . . . ,p fc+1 where npi — rc, = the number of 
observations at |,. The value of 5(1,) is then 

S&) = -4- (2.12) 

V n i 

The conditions (2.1 1) then reduce to 

ofc, + -£- ] < dfi i = 0, 1, . . . ,4+ 1 (2.13) 

v n t 

The values of n^ n u . . . ,n k+l (In,; = n) and d which give equality in (2. 13) should give a minimal value for d 
for fixed n. Solving for n. in (2. 13) we get 

* = [ , C2 ° Y (2.14) 

dy—oCi 

The solution given by eq 2.14 is a nonlinear equation in d. In cases where the sample size n is large we 
give an approximate solution which requires solving no nonlinear equations. Equation (2.13) gives 

d = -^- + -^— i = 0, 1, . . . ,4+1 (2.15) 



oc 2 



Notice that the quantities j? — = f t vary inversely withy,. 

Yi V n t 

As an approximate solution we can min max / by setting 



Hi - 



y> *(-) 2 (2.i6) 



Since by eq (2.15) d ^ — r- — this approximate solution increases the "optimum" width d by no more 

jt+i 

than c 2 o I (±.)*/y/H 

For the special form of response considered in this paper it is possible to modify the Scheffe solution to 
shrink d further. Since the resulting solution is complicated and the design can be found only by using a 
computer, this solution is presented in the appendix. 
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3. A way to implement the procedure when the knots are estimated 



It is often the case that the exact knot locations are not known precisely. However, in many situations a 
priori information is available about knot locations. In other situations there is data available so that the 
knot locations can be estimated. 

In both of the above situations the Scheffe theory does not exactly apply. However, it is reasonable to 
expect that if the estimates of the knots are good then the Scheffe theory is a reasonable description of the 
actual uncertainties involved. The exact theory for these situations is difficult and beyond the scope of this 
paper. Below we give a procedure that can be used to implement the designs derived in this paper, when the 
knots are estimated. 

The preceeding development assumes that o, £,, /?,, and y, are all known. Since in practice this is usually 
not the case, we suggest the following design plan, which incorporates a two phase procedure for collecting 
observations. 

The general design plan would then be as follows: 

(1) Take a preliminary set of m observations (these may be spaced equally or wherever they appear to give 

a good overall robust design). 

(2) With the observations from (1), estimate |„ /?,, y, and o, insert these in (2. 14) or (2.16) and solve for d so 

that Zrc, = m x > /n . 

(3) The values for n { in (2) are roughly the numbers of observations that are required at £,. Recall we already 

have m observations from (1). The remaining mr m in (2) should then be chosen to make the com- 
bined set have roughly n { observations near, and uniformly distributed about, £,. This scattering about 
4, will help to reduce potential bias caused by imperfect estimation of the knots. 

(4) Repeat steps (1), (2) and (3) if more stages are used to m 2 observations, etc. 

(5) To help implement step (3) we may proceed as follows. For a given set of observations assign each obser- 

vation to the nearest |, value. The new set of m x —m observations are then chosen as in (3). The com- 
bined set rrit are then redistributed by taking those supposedly at|, to be roughly uniform from^,-! + 
|,)/2to(|, + | m )/2. 

The remainder of this section consists of a proof of the convexity of o[c t + c 2 S(v)]. Consider the function 
S\v) defined in (2.3) and take vzh = (£,, | l+ i) for a fixed i. The only basis functions Nj[v) which are nonzero 
on /, and N^v) and N i+1 (v). Therefore 

ns 2 (v) = aN, 2 (v) + 26 N £ (v)N i+l (v) + c Nf +l (v) (3.1) 

where a, b and c are elements of the inverse of the matrix Mip) defined in (2.4). The matrix M(pi) is tridiago- 
nal, i.e. has nonzero elements m u only for \i=j\ < 1, and has only nonnegative elements. The elements a and 
c in the inverse can be seen to be positive while b is nonpositive. Equation (2.15) and some simple algebra 
then shows that (2.15) is a quadratic in v with a minimum on the interior. The convexity of o[c x + c 2 S(v)] is 
then equivalent to the convexity of 

gdO = c + d(i + (^yy" 

where c > 0, d ^ 0, toth andr ^ 0. However it can readily be shown thatg*(v) ^ 0. 

4. An illustration 

To illustrate the effect of the procedure we considered a tank which was carefully studied. The tank 
volume was approximately 13,500 L. The equation of the pressure = U versus the volume = v was thought 
to be 
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-271 + 2340 v + ll(i;—1.9K 
-27(v-4.6) + - 146(i;-5.5) + 
-12(t;-5.9) + -16(t;-6.5) + 
+ 2.2(v-7.2) + + 20(v-10.0) + 
-20(v-10.3) + + 2.5(v-12.5) + 



(4.1) 



The equation was assumed to be valid over the range | — 1700 L to | 10 = 13,500 L. Throughout the exam- 
ple we assume that the knot position, the slopes and the standard error o are known. In this case the effect of 
the design in the simplest case will be isolated. Some discussion of this will be given in the next section. The 
knots and slopes are then taken as in the above equation. The standard error o was taken to be o = 1.6 Pa. 
Using Scheffe, we then chose c x = 2 and c 2 = 5 in the equation [see (2.14)] 



* = t- 



c 2 o 



dy—aci 



Y 



(4.2) 



Using a trial and error method we solved (3.2) for d so that I! /i, = 86. (Three runs of measurements of 30 
observations each were to be taken, however 86 was more convenient than 90). The solution for d was d — 
1.03 L and the corresponding rc, values were 6.2, 6.2, 6.4, 8.2, 8.6, 8.7, 8.7, 8.4, 8.7, 8.7, 8.7. Since these were 
roughly equal it was proposed that the number of observations around each knot should be taken to be 
equal. The following two designs were then compared. 

Design 1. Take n = 86 observations equally spaced over the entire range | = 1700 to|i = 13,500. 
Design 2. Take n = 30 observations equally spaced over £ to | 10 . The remaining are chosen to make an 
approximately equal number around each knot (see part [5] in the design). 

Using the above two designs, we simulated using normally distributed error on the pressure readings with 
standard deviation 1.6 pascals for the corresponding volume readings, using eq (4.1) and ran a calibration 
chart for each design. The graphs in figure 2 show the width of the confidence intervals versus the volume 
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for the two designs. It seems clear that in the area where the knots are concentrated, between 40000 and 
8000 L, the proposed design is superior to the equally spaced design, and that the maximum for design 2 is 
significantly less than the maximum for design 1. 

5. Discussion 

The Scheffe calibration procedure involves the consideration of two bands m{v) ± o[c x + c 2 S(v)] around 
the curve u = m(v). The bands are used in a simple inversion process, where for a given reading U = w we 
solve u = m(v) for v and find an interval I(u) of possible v-values. It is required that the lengths of these 
intervals be short. A simple procedure is proposed whereby the two bands are bounded by "parallel" bands 
of uniform horizontal width and an attempt was made at minimizing the width of the outer bands. It would 
seem that a direct minimization of max I(u) would present considerable difficulty. The procedure pro- 
posed certainly needs further testing in that estimated values of | ( , /?, and o should be used in step 2 of the 
procedure. It seems, however, that the procedure will generally give lower overall width than the evenly 
spread design. 

The Scheffe'' procedure gives bands, which in our case are "square root parabolic" on segments between 
knots. In that procedure, the form of the regression function is assumed known and probabilistic statements 
are made concerning statements that the true volume i>, corresponding to a given reading U = u, is con- 
tained in an interval I(u). In our situation the unknown knots |, prevent us from assuming that our regres- 
sion functions g t (v) are known. Any attempt at making exact confidence statements is then in question. 

Since blueprints of the tanks are usually available, and n is often large, the uncertainties added by 
estimating the |, with the aid of blueprints is often not serious. However techniques such as cross validation 
(the use of half of the observations to estimate the knots and the other half to test whether these knots are 
properly located) ought to be used to assure that this is the case. 

6. Appendix 

In this appendix we take advantage of the form of the unknown function m to produce bounds which are 
shorter than Scheffe' s when o 2 is based on a large number of degrees of freedom. In order to facilitate the 
explanation we assume that o is known. 

For the curve in figure 1 the implicit bounds on m are 

- oo < ^|) - m(U <h - do 

- hi-x + c x o < mi£) - mfei) < hi - c t o i= 1, . . . ,k (A.l) 
h k+1 + c x o < mfe M ) - m(|* +1 ) < °°. 

For the Scheffe formulation of uncertainty for the calibration to hold, m must lie in the region described 
by (A.l) with probability (1-6). 

Recall that/i, = ds h and that S&) = l/v^rc,. 
In addition define 

a, = (dyt - c x o)\fni, f=0, . . . ,&+ 1 

Then the probability of m lying in the region described by (A.l) is greater than or equal to 

-^-r-(k+2)/2Uf'e- 2 Pdx) 

(A.2) 
(jV7 2 <f*)(/ e-' 1 ' 2 ) 

The correct theoretical choices of rc, £=0, ...,&+ 1 can be obtained by setting expression (A.2) equal to 
(1—6) and solving for the design which gives minimum width as in (2. 14). 
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For many purposes this solution might be difficult to obtain. Therefore it may be useful to approximate 

fe-'tpdxbyd-^—) 

Thus we choose B„ and d as small as possible so as to produce equality of 

. p -ai*f2 p ~ai2/2 p -a k + l 2 /2 

n(l -^— )(1 -*—=-. )(1 —^ ) (A.3) 

1 OiWn Otvn a k ^y n 

withl-d. 

Since the equation to be solved is nonlinear we produce an algorithm which should converge to an 
acceptable approximate solution. 

The Algorithm. 

log (A.3) if each term in the product is > 

— °° otherwise 



Define r/ . . 

f(d,n) = 



Step 1. Start with the solution (2.16) denoted by i^ and 



OCi oc 2 



mm y, fiSfni 

Step 2. Produce equality oif(d,lQ with log(l-d) by changing d if necessary, call it d . 

Step 3. Maximize/(J ,n) with respect to n and denote by n, the maximizing n, i.e. Sup/(J ,n) = f(d oy nX 
where the supremum is taken over the region n > 0, X/i, = n. 

Step 4. Iterate steps 2 and 3 generating d t and n, at the ith iteration until no significant improvement is 
possible. Round to the nearest integers ^ 1 to obtain a reasonable choice of n t . 

Proof of Convergence [for non-rounded solution]: 

Since /(d,n) is monotonically increasing in d it is clear that the sequence of rf's {d,}, obtained from the 
algorithm are monotonically decreasing. This implies that the sequence {d,} has a limit which is denoted by 
d*. 

Since the set of nonnegative rc/s such that Z n, > = n is compact it follows that every subsequent of {n,} 

has a convergent subsequence. Let {n,}, jej, be a convergent subsequence having a limit which is denoted 
byn*. 

By construction /(c?,,n,) = 1— d. Since (d j9 ikj) -* (d*,n*) asy-*«> it follows that /(<£,•-!, n,-) -* f(d*,n*) = 
1-6. 

Let !!„ denote the maximizing value of/ at d*. It follows that f(d *,!!„) ^ f(d*,n*) = 1-d. However since 
/(dj-uiij) >f(dj- u n ) ^/(d*,^) it follows thatfid*,!^) = f(d*,n*). Now since /(d* •) is strictly concave it 
follows that i^ = n*. 

Denote the optimal solution to our problem by d opt , n 09 *. 

It is clear that 

1 -d = f(d •,!!*) = /(<**', n°"). (A.4) 

The strict monotonicity of/in d implies that 

yirf*n o "0^/W op ',^ op ') 
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with equality only if d* = d opt . From eq(A.4) we have d* = d op \ and hence n* = n opt . 
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